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Abstract 

In the framework of the nonrelativistic quantum chromodynamics factorization formahsm, we 
study the processes of ip{nS) and T{nS) decay into a lepton pair or a charm pair associated with 
two jets up to the next-to-leading order in velocity expansion. We present the analytic expressions 
for the differential decay rate to the invariant mass of the lepton pair or charm pair. We find that 
the ratio of the next-to-leading order short-distance coefficient to the leading order one is in the 
range from —5.5 to —12.4. The relativistic corrections are so large that they modify the leading 
order prediction significantly. Utilizing the analytic expressions, we also investigate the relativistic 
corrections in different kinematic regions and their dependence on the masses of the initial-state 
quarkonium and the final-state fermion. In addition, we study the momentum distribution of D*^ 
in the process T(IS') — t- ccgg — > D*'^X. 

PACS numbers: 12.38.-t, 12.38. Bx, 13.20. Gd 



1 



I. INTRODUCTION 



Heavy quarkonium decay phenomena have been extensively studied both in theory and in 
experiment, from which one gains insight into both the structure of the heavy quarkonium 
and quantum chromodynamics (QCD) interactions. The predominant annihilation decay 
modes of the S'-wave spin-triplet heavy quarkonium are those hadronic decays, radiative 
decays, and leptonic decays. With abundant data of the S-wave spin-triplet heavy quarko- 
nium decays accumulated in experiments, higher order decay processes are also interesting 
to investigate. Among them, two types of processes are particularly interesting. One type 
is that the S'-wave spin-triplet charmonium and bottomonium semi-inclusive decay into a 
leptonic pair and light hadrons. The other one is the S'-wave spin-triplet bottomonium 
semi- inclusive decays into a charm meson pair and light hadrons. 

In experiment, charm production via T(IS') was studied first by the ARGUS Collabora- 
tion [1] and recently by the BABAR Collaboration [2] as well as by the CLEO Collabora- 
tion [3]. BABAR' s results [2] provided evidence for an excess of D*^ production over the 
expected rate from the virtual photon annihilation process T(1S) — ?■ 7* — ?■ cc — )■ D*^X. 
With a number of ip{nS) events accumulated at the Beijing Electron Positron Collider 
(BEPCII) [4] and T{nS) events accumulated at B factories [5], the S'-wave spin-triplet char- 
monium and bottomonium semi-inclusive decay into a lepton (charm) pair and light hadrons 
are expected to be measured well. 

In comparison with experimental data, it is necessary to theoretically study those pro- 
cesses precisely. The decay rate of these processes can be analyzed in the framework of 
nonrelativistic QCD (NRQCD) factorization formalism [6]. According to it, the decay rates 
are expressed as a sum of products of short- distance coefficients and NRQCD matrix ele- 
ments. The short-distance coefficients can be expanded as perturbation series in coupling 
constant at the scale of the heavy quark mass. The long-distance matrix elements can 
be expressed in a definite way with the typical relative velocity v of the heavy quark in the 
quarkonium state. 

The decay rate of the semi-inclusive leptonic decay process ipiT) l^l^gg was first 
studied by J. P. Leveille and D. M. Scott in the color-singlet model [7]. The polar and 
azimuthal angular distributions of the lepton pair in this process were also studied in Refs. [8, 
9]. The semi- inclusive charm decay process T — ccgg was first researched in Refs. [10, 11], 
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and the invariant mass distribution of cc has been studied in Ref. [12]. The inclusive charm 
production in T{nS) decay was calculated in Ref. [13]. Bigi and Nussinov have taken into 
account the contribution of T — )■ ccg [14]. The exclusive double charmonium production 
from T decay was calculated by Jia [15]. The authors of Ref. [16] also considered the T 
decay to two charm jets by including the color-octet contribution. Cheung, Keung, and 
Yuan calculated the color-octet J/ip production in the T decay [17]. 

According to the NRQCD factorization formula, only the leading order (LO) contributions 
are considered for the processes ip{T) {cc)gg. In the next-to-leading order (NLO), 

the decay rate receives relativistic corrections, whose long-distance matrix elements are 
suppressed by f ^ compared with the LO contribution. Notice that the relativistic corrections 
to the decay rates in the processes J/tp ^ '-/gg and J/tp ^ ggg are extremely large and 
significant [18]. One may expect that the decay rates of the processes ip{T) — l^l^gg 
and T — ccgg also receive considerable contributions from the relativistic corrections since 
those processes possess similar Feynman diagrams. However, until now, a thorough analysis 
including the contributions of the NLO NRQCD matrix elements is lacking. In this paper, we 
analyze the decay rate for ip{T) — )■ l^l~{cc)gg up to the NLO in the relativistic expansion 
in the framework of the NRQCD factorization formula. We calculate the short-distance 
coefficients of both the LO and the NLO NRQCD matrix elements at the tree level and 
present the analytic expressions for the distribution of the invariant mass of the lepton pair 
or the charm pair. With these expressions, we are able to study the relativistic corrections in 
different kinematic regions, and then provide theoretical discussions. We also investigate the 
momentum distribution of the charm quark in our work. With convolution of a charm quark 
fragmenting into a charmed hadron, we are able to predict the momentum distribution of the 
charmed hadron. Since the treatments for inclusive lepton pair production and charm pair 
production are quite similar, we concentrate on dealing with the process of H(^Si) l'^l~gg. 
The decay rate of T — )■ ccgg is readily obtained by multiplying a color factor and substituting 
the electromagnetic coupling constant a and the lepton mass mi into the strong coupling 
constant ag and the charm mass rric, respectively. 

The remainder of this paper is organized as follows. In Sec. II, we present the NRQCD 
factorization formula for the differential decay rate of the process H{^Si) gg up to the 

NLO in V. In Sec. Ill, given the notations and kinematic variables used in our calculation, we 
present the formulas for the differential decay rate as well as the total decay rate. Section IV 
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is devoted to determining the short-distance coefficients corresponding to the LO and the 
NLO NRQCD matrix elements. In Sec. V, we present the numerical results and provide 
discussions. A summary is given in Sec. VI. 



II. NRQCD FACTORIZATION FORMULA FOR QUARKONIUM DECAY PRO- 
CESS H{^Si) l+l~gg 

According to the NRQCD factorization formula, up to relative order f^, the differential 
decay rate for a quarkonium H decay into a lepton pair and light hadrons can be expressed 
as [6] 

dV[H{^S,) ^ in- +X] = {H\0,CS,)\H) + ^^^{H\V,CS,)\H) , (1) 

where m signifies the mass of a heavy quark in {H\Oi{^Si)\H) and {H\Vi{^ Si)\H) are 
the NRQCD matrix elements, and ^(^5*1) and Gi^Si) are the corresponding short- distance 
coefficients, respectively. Here H can be either charmonium or bottomonium. The four- 
fermion operators Oii^Si) and Vii^Si) are defined as 

0,{^S,) = ij^ax-X^cTij, (2a) 
1 



(2b) 



where and x Pauli spinor fields for annihilating a heavy quark, and creating a heavy 
antiquark, respectively; a* denotes the Pauli matrix; and D is the spatial part of the antisym- 
metrical covariant derivative: ip'^Dx = i^'^Dx — {D'lpVx- The subscript 1 on the NRQCD 
operator indicates that it is a color-singlet operator. According to the velocity-scaling rules 
given in Ref. [6], the matrix element of the operator Oi{^Si) in the ^Si state is of order 
while that of the operator Vi{'^Si) is of order v^. The latter one is suppressed by v"^, which 
represents the NLO relativistic corrections to the inclusive H decay. 

The vacuum-saturation approximation [6] can be used to simplify the decay matrix ele- 
ments in Eq. (2). They read 

{H\OrCSr)\H) = |(0|xV ■ e*m)\' ^ {Oi)h, (3a) 
{H\V,{'S,)\H) = Re[(iJ|^V-ex|0)(0|xV.e*(-fB)V|^^)]. (3b) 
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This approximation is valid up to corrections of relative order v^. For convenience, we 
introduce a dimensionless ratio of the vacuum matrix elements in Eq. (3) for later use [19, 20]: 



^ m^Olx^a ■ e*^\H) ' 

This quantity characterizes the typical size of relativistic corrections for H. 

Equation (1) implies that, to predict the decay rate, one needs to determine both the 
short-distance coefficients and the NRQCD matrix elements. The NRQCD matrix ele- 
ments have been extensively studied by means of lattice QCD [21], the nonrelativistic quark 
model [22], and fitting the experimental data [23, 24]. Therefore, once we determine the 
short-distance coefficients F{^Si) and G(^Si), with those values of matrix elements, we may 
calculate the differential decay rate in Eq. (1). To determine the ^(^5*1) and G{^Si) at the 
tree level, we apply the factorization formula to the process of an on-shell QQ pair near the 
threshold in a spin-triplet and color-singlet state decaying to gg: 



dV[QQ^('S^)^l+rgg] = '^^^^ {QQi('Sr)\Or{^S^)\QQ^{^S^)) 



Notice that the factorization formula (5) takes a similar form to (1) except that the hadron 
state is substituted into the on-shell free quark pair state with the same quantum number 
as the hadron. The decay rate in (5) can be calculated both in the QCD perturbation 
theory and in the NRQCD factorization formula. By matching both sides, the short- distance 
coefficients can then be determined. 



III. KINEMATICS AND FORMULAS FOR THE DECAY RATE 



A. Kinematics and definitions 



In this section, we define notations for the kinematics involved in our work. We take 
Pi and •p2 to be the momenta of the incoming heavy quark Q and heavy antiquark Q, 
respectively, which are on their mass shells: vX = V2 — ■ They are expressed as linear 
combinations of the total momentum P and half of their relative momentum q: 

pi = P/2 + q, (6a) 
P2 = P/2 - q. (6b) 
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In the center of mass frame of the quarkonium, the momenta are given by 

P = (2i?,0), (7a) 
q = (0,q), (7b) 

where the orthogonal relation P ■ g = is satisfied. 

We also assign ki, k2 to be the momenta of the two final-state gluons, and h, I2 to be the 
momenta of the produced lepton pair. Therefore, the momentum Q of the virtual photon 
yields to Q = h + I2. These momenta satisfy 

kl = kl = 0, (8a) 
ll = ll = ml (8b) 

where mi denotes the mass of the lepton. 

For convenience, we introduce a set of dimensionless variables 

2ki-P 2k2-P 2Q-P 

^1 = — ^2 = — X3 = — z = — , (9a) 

4m? |Zi I \L I I r , , , 

^- yi = Trr^ = —\-. ' (9b) 



P"^ ' \h\max mi]l 1-r'' 

where yi represents the momentum fraction for the lepton and |Zi|max denotes the maxi- 
mum of the lepton momentum in the quarkonium center of mass frame. In the following 
subsection, we will show that all the involved Lorentz invariant kinematic quantities can be 
rewritten in terms of these new variables. 



B. The formulas for the decay rate 

1. Differential decay rate of the invariant mass of the lepton pair 

For the decay process H(^Si){P) — t- {l2)g{ki)g{k2), it involves a four-body phase 

space integral, which can be expressed as 

/* - IiSkiSkiS^.iSk^''^'''''' -k.-k.~i.- (10) 

Since there is no divergence emerging in our calculation, dimensions of the space-time are set 
to 4 in (10). In order to compute the invariant mass distribution of the lepton pair, we de- 
compose the four-body phase space integral (10) into the product of a two-body phase space 
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integral for the lepton pair and a three-body one by inserting the following two identities: 

/|^(2vr) V(g -h-k) = l, P'j ^2n6iQ' - P^z) = 1. (11) 
After integrating out the energy (5° through the delta function, we get 



dz 
2^ 



(12) 



where Jd(j)2 and Jdcps are expressed as 



'--^V (2')4;(2?2,o P-''^^W (13a) 

On the other side, the squared amplitude can be expressed as the contraction of the 
leptonic tensor L'^'^ and hadronic tensor Hn^: 

\M\^ = L^^''H^,, (14) 

where the leptonic tensor L^'^ is given by 



L'^^ = 7^Tr[(fi + mi)^^{]/2 - m^Tl- (15) 



Q 

It follows after integrating over the phase space momenta that 



L^'' ^ jd^p^L^'' =[-g^'' + ^ ) X L, (16) 



where the Lorentz invariant L is given by 



2a r , r ^ , . 

where a is the fine structure constant. As a result, we are able to write the decay rate as 

r4/|i/#3S.„(-«- + ^), (18) 

where the factor | accounts for the indistinguishability of the two gluons in the final states. 
It is not hard to find that the second term in the parentheses of (18) does not contribute 
due to the current conservation. 

The three-body phase space integral fd(j)3 can generically be expressed as the integral of 
two dimensionless variables xi and X2: 

p2 r 

^3 = TOO ^ dxidx2. (19) 



Up to now, we have reduced the four-body phase space integral (10) into the integration 
over three variables: 2;, xi, and X2- The corresponding boundaries for these variables are 
given by 

^ ~ ~ ^ >X2>l-xi-z, 1 - 2 > xi > 0, l>z>r. (20) 

1 — Xi 

To simplify further the calculation, we make the variable transformation: 

xi = (1 — z)x, (21a) 
jl - z)il - x)[l - {1 - z)xy] 
- l-(l-.)x • 

After this transformation, the area of the integration is significantly simplified as 

1 > X > 0, 1 > y > 0. (22) 
Now, the expression of the decay rate reduces to 

From (23), we notice the principal task is to analyze the subprocess H(^Si) l*gg (cor- 
responding to the contribution from the hadron part H^^g^'^). In the next section, we will 
use Eq. (23) to evaluate the total decay rate as well as the differential decay rate over the 
invariant mass of the lepton pair, equivalently, the dimensionless variable z. 



2. Momentum distributions of the charm quark and the charmed hadron 

In this section, we first derive the formulas to calculate the momentum distribution of 
the charm quark in the decay process T — g*gg — ccgg. The momentum distribution of a 
charmed hadron h is then obtained by convolving it with a fragmentation function, which 
describes a charm quark fragmentation into the meson h. 

As introduced in Sec. IIIB 1, we decompose the phase space integration into two parts 
by inserting the identities (11). Since we want to observe the momentum distribution of the 
charm quark, we can integrate out the momenta of the two final-state gluons. To this end, 
we introduce a tensor T^^'^ which depends only on the momenta P and Q as 

= (-^^^ + + ^.iP' - Q'^)iP'' - (24) 
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where Hi, H2 are Lorentz invariant form factors. In the last step of (24), we have apphed 
the Lorentz covariance and current conservation. By contracting g'^'^ and pi^^P^ separately 
in (24), we are able to obtain the expressions of these two form factors. We notice that Hi 
and H2 are independent on the momenta of the two final fermions. To obtain the decay 
rate, we need to include the charm quark pair part as well as the remaining phase space. 
Contracting with the leptonic tensor,^ we readily obtain 



T - T-Z,. = i St^oStA ^-f^P - Q - A:, - A:,) X 5: X T 



(27r)32A;?(27r)32A;0' ' ^ ^ ' " m 



27rra { 



!^{r + 2z)Hi- H2 (1 - r)yl + r + z - x^^J {I - r)yl + r | 



(25) 



Now, we turn to carry out the two phase space integration /(i02 and /(i03. For /c?02, we 
have 

f \li\d\li\ 



in 



/Pi+ml\Q\ 
mf f yidy 



2vrr J ^{\-r)y\^r^x\- ^z' 



(26) 



with the boundaries of yi. 



yi+>yi> \yiA. (27) 




where 



?/i± = ^^(W1--±a/1-7)- (28) 

We then deal with the phase space integral jdcj)^. Analogously, we can reduce the integral 
Jc?03 into (19). Nevertheless, since the boundaries (28) contain X3, we prefer to choose 
another set of integration variables, such as xi and x-y. 

p2 



dx^dxi. (29) 



^ Here we should replace the leptonic tensor L^" in (15) with the corresponding tensor for the charm quark 
pair; however, we still use (15) to implement the calculation and the difference will be compensated by 
multiplying a factor. 



9 



The corresponding boundaries of X3 and xi are 



1 + z>X3> 2^^, (30a) 
Xi+ > Xi > (30b) 



where 



xi± = i (^2 - X3 ± ^xl-Az^ . (31) 

In addition, as shown in (12), to get the decay rate, we should include another integration 
over z. The corresponding boundaries of z are shown in (20) to be 1 > 2; > r. 
Finally, the decay rate can be expressed as 

^ = [dzdxsdxidy^ x T, (32) 

where T is defined in (25). 

In order to get the momentum distribution, we need to change the integration order in 
(32), and to make yi be the last integral. Notice that the boundaries of yi are independent 
of Xi; we need not change the order of the integration of xi. This calculation is tedious but 
straightforward. Here we present the expression as follows: 



" = 2 w Uo / L ^ L i L "'v L 

X xT, (33) 

V(l -r)yl + r^Jxl - Az 

where the boundaries of Xi are given in (31), and 

4± = ^ ^Jil-r)yl + r ± yi ^/il - r){z - r)z^ . (34) 
In addition, the boundaries of variable z are the positive solution of the following equation: 



l + ^±)Jl-^T(l-^± 



2± 



yi. (35) 



With the formula (33), and the boundaries (31) (34) (35), we can carry out a calculation of 
the distribution of the charm quark momentum fraction yi. Now, we go further to investigate 
the charmed-hadron momentum distribution. As discussed in Ref. [24], the momentum 
distribution of a charmed hadron produced in T decay is softer than that of the charm, due 
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to the effect of hadronization. The momentum distribution of a charmed hadron h can be 
obtained by convolving the charm momentum distribution with a fragmentation function 
for the charm quark fragmentation into the h. 

The fragmentation function Dc^h{z') describes the probabihty of a charm quark with 
hght-cone momentum /° + hadronizing into a charmed hadron h with hght-cone momen- 
tum ll + \lh\ = z'{l^ + |Zi|). The fraction z' can be expressed in terms of scaled light-cone 
momentum fractions Zi for the charm and Zh for the charmed hadron, which are analogous 
to the scaled momenta yi and yt [25], where zi is 



^1 r- — /T • l-^oj 

1 + y/1 - r 

Then, the fraction z' is expressed as 

where the last factor on the right-hand side of Eq. (37) becomes unity if the difference 
between the mass of the charm quark and that of the charmed hadron can be neglected. 
With this approximation, the momentum distribution of the charmed hadron can be written 
as [24] 

dV dzh dzi dyi dV 



dVh dyh Jzi, zi dzi dyi 



^{l-r)yl + r Jy^ \^ {I - r)yl + r + 




(38) 



where Vc-^hiz') = z'Dc^h{z'), ym. represents the upper boundary for yi in (33), which equals 



a/1 — r/2 and 1 corresponding to the first term and the second term in the parentheses, and 
Zm corresponds to the value of Zi when yi takes ym in (36). 

The formulas (33) (38), and the boundaries (31) (34) (35) can be used to carry out a 
calculation of the distribution of the charmed-hadron momentum fraction yh- In Sec. V, we 
will utilize these formulas to make predictions. 



IV. MATCHING THE SHORT-DISTANCE COEFFICIENTS UP TO NLO IN v 



In this section, we determine the differential short- distance coefficients dF{^Si) and 
dG{^Si) that appeared in (1). The short-distance coefficients are then readily obtained 
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by integrating over the integration variables. Now, we describe the strategy. By employing 
the formulas derived in the previous section, we first calculate the differential decay rate 
for the process of a color-singlet spin-triplet S'-wave heavy quark pair decay into a lepton 
pair plus two gluons QQi(^Si) — )■ l^l^gg in the QCD perturbation theory, up to the NLO 
in V, and then carry out the differential decay rate of the same process in the NRQCD 
factorization formula. Finally, the short-distance coefficients dF{^Si) and dG{^Si) in (1) 
are immediately determined by identifying these two calculations. 

A. Amplitude of QQ — ;> j*gg 

As we have demonstrated in (17) (23) (25), and (33), the lepton part has been explicitly 
written out. We still have to deal with the subprocess Q{pi)Q{p2) — ^ l*{Q)g{ki)g{k2). At 
the tree level, there are 6 diagrams contributing to the amplitude as shown in Fig. 1. Given 
the momenta defined in Sec. Ill A, the amplitude of the process reads 

A{si,S2) = v{p2,S2) Tf,u{pi,si), (39) 

where u{pi, si) and v{p2, S2) are the spinors of the heavy quark and antiquark, respectively, 
and represents the products of Dirac matrices and color-space matrices. According to 
Fig. 1, the expression of reads 

T, = {-iCQegD T'T'^® ^^(^2)7^ fi{h) ^ ^ ^ \ 7^ + 5 perms, (40) 

f2- P2-m fl+ f2- 1*2 - "1 

where e,gs denote the QED and QCD coupling constant, respectively, eg denotes the elec- 
tric charge number of the heavy quark, a, ei and b, €2 represent the color indices and the 
polarization vectors of the two gluons, and fi corresponds to the Lorentz index of the virtual 
photon. 

B. Projection of spin-triplet QQ state 

The amplitude given in (39) describes the decay of the heavy quark and the antiquark 
state with the spins of the third component si and S2, respectively. To calculate the decay 
of the QQ in the spin-triplet state and color-singlet state, one needs to project the total spin 
state of the QQ pair onto the spin-triplet and color-singlet Q{pi)Q{p2) state. This can be 
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FIG. 1: The tree-level Feynman diagrams for QQi(^Si) — ?• gg. For simplicity, the crossed 
diagrams have been suppressed. 



done by introducing the projection operator 113(^1,^2) [26] expressed by 
^3{Pi,P2) = ^u{pi,Si)v{p2,S2){-,Si;-,S2\le) ® — 



Sl,S2 

^ '\ + m){f'+2E)i{i)2-m)®^, (41) 



where Ic is the unit matrix in the fundamental representation of the color SU(3) group, 
and e is the polarization vector of the spin-triplet state. The above spin-triplet projector is 
derived by assuming the nonrelativistic normalization convention for Dirac spinor. With this 
projection operator, the amplitude for a spin-triplet and color-singlet QQ pair annihilation 
decay reads 



^ ^*9g] = Tr|n3(pi,P2)T^|, (42) 
where the trace is understood to act on both Dirac and color spaces. 



C. Projection of S'-wave amplitude 

Besides projecting the QQ pair onto the spin-triplet state, to account for the contribution 
from the S'-wave orbital-angular-momentum state, one has to project further the QQ state 
onto the S'-wave state. It can be done by averaging the amplitude over all directions 
of the relative momentum q in the QQ rest frame. 

The amplitude can be expanded in terms of the powers of and the series can be 
truncated to the desired order. Since here we are only interested in the NLO relativistic cor- 
rections, we may do it by expanding the spin-triplet amplitude A'^^^^ in through quadratic 
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order, then making the following replacement [27]: 



?V ^ yn'''^(P), (43) 



where 



pu 

U'^'^iP) ^ -g^^ + ^^. (44) 

D. The decay rate of QQi(^Si) — )• l^l~gg up to relative order v'^ 

Since the calculation for the momentum distributions of the charm quark and charmed 
hadron are similar to that of the invariant mass distribution of the lepton pair, in the 
following subsections, we merely demonstrate the latter. 

We now proceed to compute the lepton pair invariant mass distribution for the process 
QQii^Si) l^l'gg at the LO and the NLO in w, based on the techniques described in 
Sec. IV C. 

We first expand the amplitude given in (42) in terms of q up to quadratic order, then 
apply (43) to extract the S'-wave part 

A, = Af+Af^^+0{q% (45) 

The hadronic tensor in (23) is then given by squaring the amplitude A^, averaging over 
the polarizations of the initial state, and summing over the polarizations of the two gluons: 

H,. = \Y.'^X- (46) 

pol 

Substituting it into (23), the decay rate is expressed as 

4^2 , , il- zY{l-x)x ^ 1^^.^, 



2(4vr; 



, , {i-zy{i-x)x ^ 

/ dzdxdy- —— x L x - > A A 

J ^ l-(l-z)x 3 ^ 

pol 



2m2 1 /", , , (1- zf(l-x)x , 
'dzdxdy- x L 



3 (47r)4 J ^ z)x 



y + ( A^'^'A^'^* + 2Re[A^'^'A^'^*]] x 4 + 0{q') 



(47) 



X 

pol 

In the calculation, we employ the mathematica package Feyncalc [28] to implement the 
arithmetic of Dirac trace and Lorentz contraction. The resultant distribution of the invariant 
mass of the lepton pair reads 

dV Aa^aleo r r , / „ _ „ , q^ 



1-7(1 + 77-) /o(^) + /2(^)f^ , (48) 



dz 27nm'^ \l z 2z \ m? 
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where the analytic expressions for fo{z) and f2{z) are given by 
4 



{2z^ - - Viz + 8) tan"^ 



1 - z 



z 



+ 2^2(1 -^)(42^-9z + 



X tan"^ ( a/^^ ) - 9(1 - z){z^ -2z + 2) + z{hz^ - Uz + 3) logz 



(49) 



and 



9z{l - z)3 
X (17^2 -51Z + 31) 



3(4z^ - 8z^ - 57z^ + 96z - 38) tan'^ 



1 - z 



6Jz(l-z) 



tan 



;i - z){61z^ - 192^2 + 386z - 19^ 



+ 2z{z^ - 55^2 + 432 - 13) \ogz\. 



(50) 



Notice that, when extracting the relativistic corrections, we do not expand r in terms of 



E = ^J m? + in (48). Actually, from the expression of (48), we find the differential decay 
rate is sensitive to the value of r in the region of 2; — ?■ r. Moreover, the decay rate develops a 
strong dependence on r from this region, i.e., J^dz^ oc logr. In our numerical calculation, 
we will choose r = Am^/P^ = Amf/m^fj, where rriH is the mass of the initial quarkonium. 
Since the quarkonium mass is well measured, this choice may also reduce the uncertainties 
from the input parameters. 

For the same reasons, we will make the choice of r = Amj^/m'jj in (33) when evaluating 
the momentum distributions for the charm quark and the charmed hadron. 



E. The short-distance coefficients dF(^Si) and dG(^Si) 

To determine the short- distance coefficients, we need to calculate the parton level process 
QQii^Si) —7- l^l^gg in the NRQCD factorization formula. The involved matrix elements 
are easily obtained by perturbative NRQCD: 

{QQi('Si)\OiCS,)\QQ,CS,)) = 2N,, (51a) 
{QQiCSi)\ViCS,)\QQ,CS,)) = 2N,q\ (51b) 

where the state of the heavy quark pair is normalized nonrelativistically, and the factor 2Nc 
accounts for the spin and color normalization. 
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FIG. 2: Distribution of the scaled variable z for the invariant mass of the lepton pair. We use -F, G 
to signify the short-distance coefficients F(^S\) and G(^S\). 

Substituting (51) into (5), we can write down the corresponding differential decay rate 
in the NRQCD factorization formula: 

>o.f..,^.r.,^§(q:^.^s-«'^))- (-) 

Matching the QCD side and the NRQCD side by equating (48) with (52), one determines 
the short- distance coefficients '^^K and ^i^iLj^i); 

dz dz 

-/o(2)a/1 - - (1 + 7^), (53a) 



dz Slvr \ z 2z 



Employing (53), we are able to provide the following discussions. It is instructive to look 
at the ratio 

dG{'S,) , dF{^S{) Uz) 

'^'^ = -^t^'^t^- = -wr ^^^^ 

which solely depends on variable z. This ratio characterizes the importance of the NLO 
relativistic corrections compared to the LO contribution. To visualize the relation, we plot 
the ratio t{z) over the variable z in Fig. 2. From this figure, we see that the ratio t{z) is 
negative in the physical region with the variable z ranging from to 1. We also notice that 
the magnitude of t{z) rises rapidly with the increase of z. 
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In addition, we go further to analyze the two hmits of the ratio t{z). In the hmit of 
z — )■ 0, there is 

132 — IQvr^ 

which agrees with the ratio of the short- distance coefficient of the NLO relativistic corrections 
and that of the LO for the processes tp — )■ '-/gg and — ?■ ggg, as expected. In the limit of 
2 — )■ 1, it follows from Eqs. (49) and (50) that fo{z) — )■ 0, and f2{z) — )• const. As a 
consequence, 

limt(z) = -^ + ^ + 0{l-z). (56) 

From (56), we see that the ratio t{z) goes to infinity in the limit of 2 — 1, which is the 
result of a vanishing fo{z) in that limit. In fact, we can see that foiz) vanishes in the limit 
of z — 1 from amplitude. When the momenta of two real gluons are soft, the amplitude of 
J /ip ^ J* 99 can be separated into 

J, I ^ * ^ 2 /Pi ■ eiPi ■ £2 , P2- eiP2 -62 pi- eip2 ■e2+P2- eipi ■ 62 
A[J/^p ^1 gg) = g.,' 



Pi ■ kipi ■k2 P2- kiP2 ■ h Pi- kip2 ■ k2 

AiJ/^^Y), (57) 



^1112 



where and Oj indicate the polarization vector and color index of the i gluon. At LO in v, 
there is pi = P2 = and therefore A{J/ip — j- Y99) vanishes. Consequently, fo{z) vanishes 
in 2; — )■ 1. 

Figure 2 and Eq. (55) combine to indicate that the NLO relativistic corrections in this 
process are not only large but increase rapidly with the rise of the virtuality of the inter- 
mediate photon. One may doubt the convergence of the expansion series in v. In Ref. [26], 
the authors calculated the relativistic corrections to the decay rate of T — ?■ ggg up to v^. 
Their results indicate the relativistic corrections from the color-singlet matrix elements are 
convergent. Since the Feynman graphs are quite similar, we expect the relativistic expansion 
will be convergent in the process of T — )■ l^l^gg. 

The short-distance coefficients G(^Si) and F{^Si) can be readily obtained by integrating 
out the variable z. Finally, substituting the short-distance coefficients given in Eqs. (53) 
into Eq. (1), we present the differential decay rates in the NRQCD factorization formula for 
the process H(^Si) — t- l^l~gg: 
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TABLE I: Numerical values for the parameters of different initial-state particles: the mass rriH, 
strong coupling constant as{mi{ /2), the NRQCD matrix elements {Oi)h and the value of {v'^)h- 





ruH {GeY) 




(Oi)j^(GeV3) 




J/i; 


3.097 


0.334 


0.440 


0.225 




3.686 


0.300 


0.274 


0.633 


T(15) 


9.460 


0.215 


3.07 


0.057 


T{2S) 


10.023 


0.211 


1.62 


0.179 


T{3S) 


10.355 


0.210 


1.28 


0.251 



Tz = + (^^^ 

where the matrix element {v'^)h is previously defined in (4). The decay rate is correspond- 
ingly achieved by integrating out the variable z. 

The differential short- distance coefficients as well as the decay rate for the lepton pair 
production can be easily extended to the process T{nS) — )■ ccgg, where the charm pair 
is produced through one virtual gluon instead of the virtual photon. One can get them 
by multiplying a color factor 5/24, and substituting mi and CqO^ into niD and on the 
right-hand side of (53) and (58).^ 

V. NUMERICAL RESULTS AND DISCUSSIONS 

In this section, we first numerically evaluate the total decay rate and the short-distance 
coefficients for various processes, and then discuss the momentum distribution related to 
the charmed meson D*^ in the process T(IS') — ?■ ccgg — )■ D*^X. 

A. Decay rate and the short-distance coefficients 

In this subsection, we employ the obtained differential short-distance coefficients (53) 
and the decay rate (58) to make numerical predictions for the decay rate of the processes 

^ In calculating the decay rate of the process TinS) — > ccgg, we take the mass of the charm quark to be 
that of the D meson nic — mn in order to compare with the measurement of the experiment [13]. 
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TABLE II: The ratio r, theoretical predictions for the decay rate, the ratio between the NLO rate 
and LO rate, and the ratio between the short-distance coefficients. 







r 




r(o)(keV) 


r(2)(keV) 




r(2) /r(o) 


Gi(^Si)/Fi{^Si) 




e+e gg 


1.08 X 10" 


-7 


4.73 X 10"^ 


-5.91 X IQ- 


-1 


-125% 


-5.56 


J/iJ ^ 


fi~^H~gg 


4.69 X 10" 


-3 


1.08 X 10-1 


-1.57 X 10" 


-1 


-145% 


-6.49 


^(25) ■ 


-)■ e~^e~gg 


7.66 X 10" 


-8 


2.43 X 10-1 


-8.56 X IQ- 


-1 


-352% 


-5.55 


^{2S) ■ 


-?■ i-i^fi^gg 


3.31 X 10" 


-3 


5.98 X 10-2 


-2.41 X 10" 


-1 


-403% 


-6.37 


T(1S) 


e~^e~gg 


1.16 X 10" 


-8 


3.68 X 10""2 


-1.16 X IQ- 


-2 


-31.5% 


-5.53 


T(1S) 




5.02 X 10" 


-4 


1.22 X 10-2 


-4.16 X 10" 


-3 


-34.0% 


-5.97 


T(15) 




1.41 X 10" 


^1 


1.05 X 10-3 


-7.06 X IQ- 


-4 


-67.3% 


-11.8 


T{IS) 


ccgg 


1.56 X 10" 


-1 


1.44 


-1.01 




-70.4% 


-12.4 


T(25) 


ccgg 


1.39 X 10" 


-1 


7.99 X 10-1 


-1.68 




-210% 


-11.7 


T(35) 


ccgg 


1.30 X 10" 


-1 


6.63 X 10-1 


-1.90 




-287% 


-11.4 



H{^Si) —J- {cc)gg. The corresponding discussions are also presented. 

To this end, we need to specify various input parameters, such as the couphng constants, 
the pole masses of the heavy quarks, the physical masses of various involved quarkonia and 
final-state leptons and charm quark (we choose the mass of the final-state charm quark 
to be the mass of the charmed hadron), and the values of the nonperturbative NRQCD 
matrix elements. In our calculation, we take the charm and bottom quark pole masses to 
be rric = 1.4 GeV and m;, = 4.6 GeV, respectively. The lepton masses are taken to be 
me = 0.51 X 10-3 GeV, = 0.106 GeV, = 1.777 GeV [29]. Since the final-state charm 
quark will dominantly evolve to the charmed hadron, we choose the charm quark mass to 
be the mass of the charmed hadron rriD = 1.87 GeV, which is the average masses of the 
and The fine structure constant changes slightly from the scale of charmonium to that 
of bottomonium, so we uniformly choose a = ^ for all the decay processes involved. 

The values of the quarkonium masses, coupling constants, and the NRQCD matrix ele- 
ments are listed in Table I, where scales of the coupling constants are chosen to be half of the 
corresponding decay quarkonium. In the table, the masses of the quarkonia are taken from 
Ref. [29]; we take the NRQCD matrix elements and {v'^)j/^ from Ref. [24], (O)T(nS) 

from Ref. [12], and {0)^{2S) from Ref. [30]; other values of the NRQCD matrix elements 
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{v'^)h are determined by the Gremm-Kapustin relation [19]: 



3 



{v^)h = "^"^ (59) 

where vripoie denotes the pole mass of the heavy quark, which is taken to be 1.4 GeV and 
4.6 GeV for the charm quark and bottom quark, respectively. 

With the parameters chosen above, we are able to make numerical predictions for various 
decay channels, which include the inclusive lepton decay of the charmonium and bottomo- 
nium, as well as the inclusive charm decay of the bottomonium. First, we consider the total 
decay rate. The predicted results are listed in Table II. In the table, we give the decay 
rates both in the LO and in the NLO relativistic corrections. To show the magnitude of 
the relativistic corrections, we also list two ratios. One is the ratio of the NLO and the LO 
short-distance coefficients, namely, G(^S\)I F{^S\). The other is the ratio of the NLO and 
the LO decay rates r(2)/r(°). 

From Table II, we find that all the relativistic corrections are huge and negative. This 
is especially serious for the bottomonium decay to charm pair channels. We can reach two 
conclusions from the table. First, the ratio of the NLO and the LO short- distance coefficients 
ascends with the increase of r, which is previously defined as 4m^/m^ [or 4m|)/m^ for 
T(nS') — ccgg\. Second, in the channel with small r such as T(?t,S')(^/'(?t,S')) — )■ e^e~gg, the 
ratio of the short- distance coefficients approaches to that in the process of J/ip ^ 'ygg or 
J/tp ^ ggg. This is understood from the fact that the decay rate of H{^Si) — > e~^e~gg is 
dominated by the region, where the virtual photon is nearly on-shell. 

It is also intriguing to study the r dependence of the relativistic corrections. In Fig. 3, we 
show the dependence of the ratio Gl'^Si)/ F(^Si) on r. From the figure, we see that as the 
value of r increases, the relativistic corrections will increase rapidly. Actually, this feature 
has been shown in Table II. When the mass of the final-state fermion is close to half of 
that of the initial quarkonium, the momenta of the two real gluons will become soft, and 
therefore the perturbative QCD calculation is unreliable. Therefore, only the region r < 0.5 
is plotted in Fig. 3. 



^ Since the pole masses of the charm quark and bottom quark are not determined very well, the NRQCD 
matrix element computed from the Gremm-Kapustin relation has a large uncertainty. This is especially 
serious for the bound state quarkonium, whose mass is close to 2mpoie- Therefore, in the next subsection, 
we adopt a new method to determine {v'^)h for T(IS'). 
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-10- 
£ -15- 
O -20- 

-25- 

-30- 

FIG. 3: Dependence of the ratio of the short-distance coefficients G{'^Si)/F{^Si) on r. 
B. Momentum distribution of charmed hadron D*^ 

To predict the production rate of a charmed hadron from T(IS') decay, we need to consider 
the probabihty of a charm quark hadronizing into the charmed hadron. In Ref. [31], the 
authors computed the ratio Br[c h]. In the Table 10 of Ref. [31], one can read that the 
ratio for D*'^ production is Br[c — t- D*^] = 0.220. With this value, we can readily derive 
the decay rate for D*"^ production through the process T(IS') — )■ ccgg — )■ D*~^X. 

As mentioned in the previous subsection, the NRQCD matrix element {v'^)r{is) deter- 
mined from the Gremm-Kapustin relation is sensitive to the bottom pole mass. Here we 
present another method to determine this matrix element, and then use the new value to 
predict the momentum distribution of D*^. 

In Ref. [2], the BABAR Collaboration reported their measurement Br [T(IS') — i- D*'^X] = 
(2.52±0.13(stat) ±0.15(syst))%. In addition, they derived the contribution from the virtual 
photon annihilation process to be Br [T(IS') — 7* — D*+X] = (1.52±0.20)%, and therefore 
we may expect that the difference arises from the contribution of T(IS') — > ccgg — D*~^X. 
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TABLE III: The KLP and Peterson fragmentation function and the value of the corresponding 
parameters. 





D{z') 


Nh 




KLP 




n.o 


5.6 


Peterson 




-2 0.127 


0.054 



With this assumption/ we are able to fix the value of (f^)T(i5) through the relation^ 

^(f{^S^) + G{^S^){v^)i^ X X Br[c^ = 2.52% - 1.52% = 1.00%, (60) 

i T V J 

where Fx represents the total decay rate of T(15'). By taking as Fx = 54.02 keV, we obtain 

{v'^)t{is) = —0.0781. In Ref. [32], this matrix element is also determined to be — 0.009±0.003 

through fitting the decay rate of the process T — )■ e^e~ . Though both results are negative, 

our result is much larger than theirs. We apply the formulas (33) and (38) derived in 

Sec. IIIB2 to calculate the momentum distribution of D*'^ . Prior to making the numerical 

predictions, we need to select an appropriate fragmentation function. Here we employ two 

well-known models: the Kartvelishvili-Likhoded-Petrov (KLP) fragmentation function [33], 

which was used in the analyses of charmed-hadron momentum distribution in T{nS) and Xh 

decays, and the Peterson fragmentation function [34]. The KLP and Peterson fragmentation 

functions both have a simple parametrization depending only on the light-cone momentum 

fraction z' (see Table III). The optimal values of determined by the Belle Collaboration 

are Uc = 5.6, and 0.054 for the KLP and Peterson fragmentation functions, respectively [31]. 

The normalization factor Nh is determined by the constraint dzDc^h{z) =Br[c — )■ h]. 
Taking the fragmentation probability Br[c — )> D*^] = 0.220, we are able to determine the 
normalization factors for the two fragmentation functions, which are shown in Table III. 

With the formulas (33) (38) and the fragmentation functions in Table III, we can evaluate 
the momentum distribution of D*~^ in the KLP and Peterson models, which is shown in 

^ In Ref. [16], the authors considered the contribution to the charm pair production from the eolor-octet 
NRQCD matrix element. According to the NRQCD velocity-scaling rules, this contribution belongs to the 
higher order corrections at v expansion. We are now working on corrections to the process T ccgg, 
and a thorough analysis including the contributions from the color-octet NRQCD matrix elements will be 
presented in the future. 

^ Since we use the experimental data related to the D*^ production, here we choose r = Am\i,+ /mi^^-^gy 
where m^.+ = 2.01 GcV. 
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Fig. 4. We notice that the discrepancy between the figures from the two models is small. 
This implies the momentum distribution is insensitive to the models. In addition, we find 
that the contribution from the NLO relativistic corrections is comparable with that of the 
LO, and therefore modifies the LO magnitude significantly. 




FIG. 4: The momentum distribution of the charmed hadron D*~^ for the different fragmentation 
function. The left figure is for the KLP model, and the right one is for the Peterson model. In 
the figure, the dotted line, dash-dotted line, and solid line correspond to the LO, NLO, and total 
distributions, respectively. 

VI. SUMMARY 

In this work, we compute the NLO relativistic corrections to the decay rates of the 
processes of ip{nS){T{nS)) — t- l'^l~{cc)gg in the framework of the NRQCD factorization 
formula. The differential short-distance coefficients and decay rates over the invariant mass 
of the lepton pair or the charm pair are presented analytically. The relativistic corrections 
to all the processes are significant. The magnitude of the NLO relativistic corrections even 
surpasses that of the LO contribution in most processes. Furthermore, we analyze the 
ratio of the differential short-distance coefficients. The results indicate that the relativistic 
corrections increase rapidly with rise of the invariant mass of the lepton pair or the charm 
pair. In addition, we study the r dependence of the ratio of the short-distance coefficients 
G(^Si)/F{^Si). In the limit of r — 0, our result is consistent with that of J/ip — )■ ^gg or 
J/ip ^ ggg. With the increase of r, the ratio G^^Si)/ Fi^^Si) increases rapidly. 

The momentum distributions of a free charm quark and of a charmed hadron in the 
process T(IS') — ccgg — )■ DX are studied. We also determine the NRQCD matrix element 



23 



(f^)x(is) based on the measurement of the BABAR Collaboration. Taking it as an input 
parameter, we also predict the momentum distribution of D*^ through the process T(IS') — )• 
ccgg D*+X. 
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